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A Pseudospectral Method for Optimal Control 
of Open Quantum Systems 
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In this paper, we present a unified computational method based on pseudospectral approxima- 
tions for the design of optimal pulse sequences in open quantum systems. The proposed method 
transforms the problem of optimal pulse design, which is formulated as a continuous time optimal 
control problem, to a finite dimensional constrained nonlinear programming problem. This resulting 
optimization problem can then be solved using existing numerical optimization suites. We apply the 
Legendre pseudospectral method to a series of optimal control problems on open quantum systems 
that arise in Nuclear Magnetic Resonance (NMR) spectroscopy in liquids. These problems have 
been well studied in previous literature and analytical optimal controls have been found. We find 
an excellent agreement between the maximum transfer efficiency produced by our computational 
method and the analytical expressions. Moreover, our method permits us to extend the analysis and 
address practical concerns, including smoothing discontinuous controls as well as deriving minimum 
energy controls. The method is not restricted to the systems studied in this article but is universal 
to every open quantum system whose performance is limited by dissipation. 



PACS numbers: 



I. INTRODUCTION 



The problem of relaxation is ubiquitous in all applica- 
tions involving coherent control of quantum mechanical 
phenomena. In these applications, the quantum system 
of interest interacts with its environment (open quantum 
system) and relaxes back to some equilibrium state 1]. 
This relaxation effect leads to degraded signal recovery 
and, in turn, to the loss of experimental information. Op- 
timal manipulation of open quantum systems in such a 
way as to produce desired evolutions while minimizing re- 
laxation losses has been a long-standing and challenging 
problem in the area of quantum control. 

Various methods employing optimization techniques 
and principles of optimal control have been developed 
for the design of pulse sequences that can be used to ma- 
nipulate quantum systems in an optimal manner. How- 
ever, the large majority of them is limited to deal with 
closed quantum systems [1, H, 0, [E IE S Hi- Recently, 
relaxation-optimized pulse sequences that maximize the 
performance of open quantum systems have emerged. In 
some simple cases, for example maximizing polarization 
transfer between a pair of coupled spins in the presence of 
relaxation, the optimal pulses have been derived analyt- 
ically using optimal control theory [l^, |Tl[ . For more 
general cases, gradient ascent algorithms were proposed 
to optimize pulse sequences for optimally steerin g th e dy- 
namics of coupled nuclear spins [ij, [isl, [3, [iM IH) [13 ■ 
These algorithms, while successful, rely heavily on the 
computation of an analytic expression for the system evo- 



lution propagator and gradients as well as a large number 
of discretizations over which to evolve the system. This 
results in expensive computational power, and gradient 
ascent algorithms in general inherit slow linear conver- 
gence rates (18i] . 

In this article, we present a unified computational 
method for optimal pulse sequence design based on pseu- 
dospectral approximations. This paper is organized as 
follows. In the next section, we formulate optimal control 
problems in open quantum systems and introduce pseu- 
dospectral methods. In Section IIIIl we present several 
examples to demonstrate the robustness of pseudospec- 
tral methods for optimal pulse sequence design. The sys- 
tems in the examples have been thoroughly studied in 
our previous work or by others. 



II. PSEUDOSPECTRAL METHODS FOR OPEN 
QUANTUM SYSTEMS 

For an open quantum system, the evolution of its den- 
sity matrix is not unitary. In many applications of inter- 
est, the environment can be approximated as an infinite 
thermostat the state of which never changes. Under this 
assumption, so-called the Markovian approximation, it is 
possible to write the evolution of the density matrix of 
an open system (master equation) alone in the Lindblad 
form fii 



p=-i[H{t),p\-L{p) , (ft=l). 



(1) 
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where H{t) is the system Hamiltonian that generates 
unitary evolution while the term L{p) models relaxation 
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(nonunitary dynamics). The general form of L is 

Q,/3 



(2) 



where Va^ p are operators that represent various relax- 
ation mechanisms and kap are coefficients that depend 
on the physical parameters of the system. The Hamilto- 
nian H{t) has the general form 



H{t)=Hf + Y,Mt)H^ 



(3) 



where Hf is the free evolution Hamiltonian and 
^™-^Ui(t)_ffi is the so-called control Hamiltonian. The 
later is used to manipulate the open quantum system by 
application of electromagnetic pulses Ui(t) of appropriate 
shape and frequency. 

A typical problem of controlling an (finite-dimensional) 
open quantum system has the following form: starting 
from some initial state p{0) at i = 0, find optimal pulses 
Ui{t) , < t < T , that bring the final density matrix p(T) 
a.t t — T as "close" as possible to some target opera- 
tor O. More precisely, find Ui{t) that maximize the final 
expectation value of O, {0){T) = tra.ce{p{T)0}. Us- 
ing the master equation ^ and the relation d{0)/dt — 
tracej/jO} that holds for a time-independent operator, 
we can find a system of ordinary differential equations 
that describe the time evolution of an open system for 
the desired transfer p{0) O 20] 



(4) 



where x = (xi, . . . , a;„)^ G M" is the state vector whose 
elements are expectation values of the operators par- 
ticipating in the transfer (for example usually Xn{t) = 
{0){t)), HfjHi € R"^" are square matrices correspond- 
ing to operators Hf,Hi under a fixed basis of the state 
space (Hilbert space) and t g [0, T]. This gives rise to 
an optimal control problem that starting from an initial 
state x{0) (which is related to p{0)), find the controls 
Ui{t), t g [0,r], that maximize Xn{T) = (0)(T) sub- 
ject to the system evolution equations as in Specific 
examples are given in the next section. 

Practical considerations such as power and time con- 
straints guide us to consider a more general cost function 
(the quantity that we want to maximize or minimize) 



min ip(T,x(T)) 



£{x{t),u{t))dt , 



(5) 



where u — (wi,U2, . . . ,Um)^ is the control vector, (p is 
the terminal cost depending on the final state at the ter- 
minal time t = T, and C is the running cost depending 
on the time history of the state and control variables, x 
and u. For example, if ip = —Xn(T) and £ = 0, then ^ 
is equivalent to maximizing Xn{T) as mentioned above. 



while if (p = and C = YllLi'^'i^ then ([5]) is equivalent 
to minimizing the energy of the pulses. In many cases of 
application, not only the initial state a;(0) can be speci- 
fied but also other endpoint constraints may be imposed. 
They can be expressed in a compact form as 



z{x{Q),x{T))=Q. 



(6) 



Additionally, constraints on the state and control vari- 
ables satisfied along the path of the system may be im- 
posed, such as the amplitude constraints where |'Ui(t)| < 
M , for all t S [0, T], where M is the maximum amplitude 
of the pulses. Such constraints can be expressed as 



g{x{t),u(t))<Q. 



(7) 



This class of optimal control problems of bilinear systems 
(linear in both state and control) described in (j4]) are in 
general analytically intractable. However, they can be 
efficiently solved by pseudospectral methods. 

Spectral methods involve the expansion of functions 
in terms of orthogonal polynomial basis functions on the 
domain [—1,1]. Using such a basis leads to spectral accu- 
racy, namely, the k^^ coefficient of the expansion decays 
faster than any inverse power of k |21], which is anal- 
ogous to the Fourier series expansion for periodic func- 
tions. This property of rapid decay from spectral meth- 
ods is adapted to solve optimal control problems like the 
one described above. It permits the use of relatively low 
order polynomials to approximate the control and state 
trajectory functions, u{t) and x{t). 

Since the support of the orthogonal polynomial bases 
is on the interval [—1, 1], we first transform the optimal 
control problem from the time domain t € [0, T] to r g 
[—1, 1] using the simple affine transformation. 



2t-T 
T ' 



In a redundant use of notation, we make this transition 
and reuse the same time variable t. The transformed 
optimal control problem can now be written as, 

T 

min (/7(l,a;(l)) + - / £{x{t),u{t))dt 
T r ™ 

s.t. Hf + ^u,{t)Hi X (8) 

i=l 

e(a;(-l),x(l)) = 
g{x{t),u{t)) < 0. 

Pseudospectral methods were developed to solve partial 
differential equations and recently adapted to solve opti- 
mal control problems 22, 23, 24, 25, 26]. Pseudospectral 
approximations are a spectral collocation (or interpola- 
tion) method in which the differential equation describing 
the state dynamics is enforced at specific nodes. Spectral 
collocation is motivated by the Chebyshev Equioscilla- 
tion Theorem [27| which states that the best N^^ order 
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approximating polynomial p*pf{f) to a continuous func- 
tion / on the interval [—1, 1] is an interpolating polyno- 
mial, as evaluated by the uniform norm, 



II/-Pa'(/)IIoo = mill 



\f-p\\ 



(9) 



where Pn is the space of all polynomials of degree at most 
N. Since any iV**^ order interpolating polynomial can be 
represented in terms of the Lagrange basis functions (or 
Lagrange polynomials), we use these functions to express 
the interpolating approximations of the continuous state 
and control functions, x{t) and u{t), as in the model ([8]). 
Given a grid of -I- 1 interpolation nodes within [—1,1], 
r = {^0 < ti < ■ ■ ■ < ijv}: the Lagrange polynomials 
{ik} e Pn, fc e {0, 1, . . . , N}, are constructed by 



N 



= n 



which are characterized by the /c*'' polynomial taking unit 
value at the k^^ node of the grid and zero value at all 
other nodes of the grid, i.e., ik (U ) = 5ki, where Ski is 
the Kronecker delta function [2g|. Note that Lagrange 
polynomials form an orthogonal basis with respect to the 
discrete inner product {p,q) — X]fcLo?'(*fc)'?(^fc)- 

With these tools, we can now write the N^^ order in- 
terpolating approximations of the state trajectory and 
control functions with respect to a given grid F of -|- 1 
nodes as, 




-0.5 0.5 

X 

(a) Uniform Grid 




-0.5 0.5 

(b) LGL Grid 



FIG. 1: The A'' = 16 order interpolation of the function 
f{x) = l/(16x^ + 1) based on a uniform grid (panel a) 
demonstrates the Runge phenomenon whereas the interpo- 
lation based on the LGL grid (panel b) does not. 



product, 



{f,9) = I J{t)g{t)w{t)dt, 



with a constant weight function w{t) — 1, Vt S [—1, 1], 
where f,g € L2[— 1, !]• Implementing the pseudospectral 
method with the orthogonal Legendre polynomials deter- 
mines the grid to be Legendre-Gauss nodes which are the 
roots of Ljv(i), the derivative of the iV'^ order Legendre 
polynomial. To enforce the method at the end points, we 
use the Legendre-Gauss-Lobotto (LGL) nodes, which in- 
clude to = -1 and In = 1, i.e., T^'-'^ = {tj : LN{t)\t^ = 
0, j = 1, . . . iV - 1} U{-1, !}• Then, the Lagrange poly- 
nomials ik{t) can be expressed with respect to p-^*^^ as 



x{t) « lNx{t) = Y.k=o ^kf-k{t), 



i{t) « lNu{t) Yl,k=a Wfe4(^), 



(10) 
(11) 



where Xk and Uk are not only the coefficients of the ex- 
pansions, but also the function values at the k^^ node due 
to the definition of the Lagrange polynomials [l^]. Be- 
cause these Lagrange polynomials are constructed based 
on the choice of these nodes, the approximations made 
with this basis as in pH)) and pT|) are sensitive to the 
choice of the nodes. For an arbitrary selection of nodes, 
as the order of approximation N gets large, Runge phe- 
nomenon may occur, that is, there are increasingly larger 
spurious oscillations near the endpoints of the [—1, 1] do- 
main [2!|| as shown in Fig. [T] A selection of Gauss-type 
nodes with quadratic spacing towards the endpoints sup- 
presses such oscillation between the interpolation nodes 
and greatly increases the accuracy of the approximation 
[30| . It has been shown that for a fixed > and a 
norm given by (O, Gauss- type nodes are asymptotically 
close to optimal for interpolating a continuous function 
over the domain [—1, 1] [3lj . 

In order to maintain the advantages of a spectral 
method while using collocation, we write the Lagrange 
polynomials in terms of orthogonal polynomials. We 
choose to focus on the Legendre polynomials which are 
orthogonal in L2[— 1,1] defined with a weighted inner 



{e-i)LM{t) 



N{N + l)LN{tk) t-tk 



where {tk} € F^'^^, fc = 0, 1, . . . , TV. 

From the interpolation as in pop . we have 



N 



— lNx{t) = ^ Xklkit). 



k=0 



Using special recursive identities for the derivative of 
Legendre polynomials [1^, we have at the LGL nodes 



N 



N 



dt 



lNx{tj) = ^Xkikitj) = ^DjkXk, (12) 



k=0 



fe=0 



where Djk are jk elements of the constant (TV -I- 1) x 
(A^ -1-1) differentiation matrix D defined by [S^l 



Djk - <^ 



jV(jV+l) 



jV(JV+l) 
4 



j = k = 
j = k = N 



otherwise. 



(13) 
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In addition, the integral cost functional in the opti- 
mal control problem ([S]) can be approximated by Gaus- 
sian quadrature. In particular, Legendre-Gauss-Lobotto 
quadrature is used to enforce endpoint conditions and 
defined as 



^1 N 1 

J f{t)dx = ^ f{ti)wi, J 



ir{t)dt, (14) 



which is exact for / G P2Ar-i when {U} e T^'^^ [2l|. 
Therefore, the choice of LGL nodes not only achieves 
close to optimal interpolation error by preventing increas- 
ingly spurious oscillations as N gets large but also ensures 
the accuracy of the numerical integration. 

Gompihng equations (fTO|) . (fT2|) . and (fT4|) we can con- 
vert the optimal control problem as in ^ into the fol- 
lowing finite-dimensional constrained minimization prob- 
lem by discretizing the states and controls with an inter- 
polation scheme, representing the differential equation 
through recursive definition of spectral derivatives, and 
expressing integral terms with Gaussian quadrature. 



min Lp{T, xpf) 



T 



N 



[Xi,Ui)Wi 



N 



s.t. 



^3 ' 



^^^^^ ^ 

fc=0 ' i=l 

e(xo, xn) = 0, 

5(x,,w,)<0, VjG{0,l,...,iV}, 



where Uij,i = l,...,r7i, are components of the vec- 
tor Uj denoting the value of the control function ui at 
the j'^ LGL node t,, namely, Uj = {uij, . . . ,UmjY' — 
{ui{tj), . . . , Um{tj)) ■ Solvers for this type of constrained 
minimization problem are readily available and straight- 
forward to implement. 



III. EXAMPLES FROM NUCLEAR MAGNETIC 
RESONANCE SPECTROSCOPY IN LIQUIDS 

In this section we show the robustness and efficiency 
of the pseudospectral method by applying it to a series 
of optimal control problems on open quantum systems 
that arise in NMR spectroscopy of proteins in liquids. 
These control problems were selected because analytical 
expressions for their optimal solutions have been derived 
in the literature d, [lO, [I3| , making them well suited for 
testing the performance of the pseudospectral method on 
open quantum systems. 



A. Pair of Coupled Heteronuclear Spins 

The first open quantum system from liquid state NMR 
that we consider is an isolated pair of heteronuclear spins 



1/2 (spins that belong to different nuclear species), de- 
noted as Ii (for example ^H) and I2 (for example ^'^C or 
^^N), with a scalar couphng J [13]. In a doubly rotating 
frame, which rotates with each spin at its resonance (Lar- 
mor) frequency, the free evolution Hamiltonian for this 
system is Hf = 2JIiJ2z, where /i^ = aiz/2, hz = ^22/2 
and uiz, (J2z are the Pauli spin matrices for spins Ii and 
I2 respectively. Note that this Hamiltonian is valid in 
the so-called weak coupling limit, where the resonance 
frequencies of the spins satisfy \loi — W2I ^ J and thus 
the Heisenberg coupling (Ji • I2), which is the character- 
istic indirect coupling between spins in isotropic liquids, 
can be approximated by the scalar coupling (Iizhz) [lOl- 
The most important relaxation mechanisms in NMR 
spectroscopy in liquid solutions arc due to Dipole-Dipole 
(DD) interaction and Chemical Shift Anisotropy (CSA), 
as well as their interference effects, i.e., DD-CSA cross 
correlation ^^l- We initially consider the spin system 
without cross-correlated relaxation. 



1. Spin Pair without Cross- Correlated Relaxation 

Here we consider the open quantum system with only 
DD and CSA relaxation ignoring the cross-correlated re- 
laxation. This case approximates for example the situa- 
tion for deuterated and ^^N-labeled proteins in H2O at 
moderately high magnetic fields (e.g., 10 Tesla), where 
the ^H-^^N spin pairs are isolated and CSA relaxation 
is small. Furthermore, we focus on slowly tumbling 
molecules in the so-called spin diffusion limit [3^. In 
this case, longitudinal relaxation rates (1/Ti) are negli- 
gible compared to transverse relaxation rates (I/T2) f33l|. 

For this coupled two-spin system, the free evolution of 
the density matrix p in the doubly rotating frame is given 
by the following master equation 9] 

p = - iJ[2Iizl2z,p] - kooWizhz, Wizhz,p]] 

- kcsAihz, [hz,p]] - kl'SAihz, [l2z,p]], (15) 

where J is the scalar coupling constant, koD is the DD 
relaxation rate, and fc^^^, fc^SA CSA relaxation rates 
for spins /i , /2 , respectively . These relaxation rates de- 
pend on various physical parameters of the system, like 
the gyromagnetic ratios of the spins, the internuclear dis- 
tance and the correlation time of the rotational tumbling 

One problem of interest in NMR experiments is to find 
the pulses (controls), LUxit) and ujy{t) applied in the x and 
y directions, respectively, for optimal polarization trans- 
fer Iiz I2Z from one spin to the other. This trans- 
fer is suitably done in two steps: Iiz — > 2Iizl2z Iiz- 
Since these two steps are symmetric, the optimal con- 
trols for the second step are symmetric to those of the 
first one. Thus, we only need to concentrate on the 
first step and the objective is to maximize the trans- 
fer Iiz — > 2Iizl2z- In other words, starting from the 
initial state p(0) = Iiz, we tend to maximize the final 
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expectation value of the target operator O = ilizhz, 
i.e., (2/i,/2,)(T) = trace{p(r)2Ji,/2j, where T is the 
final time of the experiment and p{T) is controlled by 
Lu,j:{t) and LOy{t). Using the master equation (|15p . we can 
find differential equations that describe the time evolu- 
tion of the expectations of the operators participating 
in the desired transfer as presented in Section |TT1 The 
corresponding equations in matrix form are 



(16) 



where xi = {h^), X2 = {hx), x^ = {2Iiyl2z), x^ = 
(2/1^/22), ^ = (fcuD + fccS/i)/"^' ^"^^ controls Ui{t) = 
LOy{t)/J and U2{t) = uJx(t)/J are the normalized (w.r.t. 
J) transverse components of the applied magnetic field. 
Note that the above system ([TBI) the bilinear form as 
shown in 

Consequently, we now arrive at an optimal control 
problem for the transfer Ii^ "^Iizhz that is to find 
ui{t) and U2{t), < t < T, such that starting from 
^(0) = (1,0,0,0)-^, X4(T) is maximized subject to the 
evolution equations (fTB|) . This problem has been solved 
analytically and the resulting analytical pulse was de- 
noted as ROPE [9] . It is shown there that the maximum 
achievable value of 0:4, i.e., the efficiency 771 of the transfer 
is given as a function of parameter ^ by 



m 



(17) 



Using the pseudospectral method presented in this 
paper, we calculated numerically optimal controls 
ui{t),U2{t) for various values of ^ in the range [0,1] to 
maximize the corresponding achievable values of x^iT). 
The method was implemented in Matlab using the third 
party KNITRO nonlinear programming solver from Ziena 
Optimization. The problem was approximated using 25 
{N = 24) nodes and with the terminal time free to vary, 
but with a maximum time of T = 10. A unified and gen- 
eral method should not use any prior knowledge, so the 
solver was given an arbitrary initial guess for the controls. 
In this case, we take ui{t) = U2{t) — 1 and T — 1. The 
termination tolerance on the cost function of the solver 
was set at 1 X 10~^. In Fig. [5] we plot the value X4(T) 
achieved by the pseudospectral method for ^ S [0, 1]. For 
comparison, we also plot the maximum efficiency given 
in (fTT)) . The excellent agreement shows the efficiency of 
the method to approximate optimal solutions. 

Another clear advantage of the pseudospectral method 
well illustrated by this problem is that the derived control 
pulses are smooth functions. Fig. 3(a) shows the discon- 



tinuities of the theoretically calculated optimal pulse am- 
plitude (xA^iW+mIW) M- Such discontinuities can be 
challenging, if not impossible, to implement in practice 
and high amplitudes can be hazardous for the experiment 
sample, equipment, and human subjects as in Magnetic 
Resonance Imaging (MRI). The pulse amplitude derived 



by the pseudospectral method, shown in Fig. |3(b)[ is 
easily implementable and maintains low values despite 
achieving transfer efficiencies within 1 x 10^^ of the the- 
oretical optimal values. The pseudospectral pulse shown 
in Fig. |3(b)| is attained from an optimization that mini- 
mizes energy while maintaining a desired efficiency trans- 
fer (i.e., the value of 2:4 (T) found without taking energy 
considerations into account). Since there are many pulses 
that achieve this near-optimal transfer this optimization 
selects the pulse with minimum energy. Therefore, not 
only is the pseudospectral pulse without discontinuities 
but it also accomplishes the transfer with 45% less energy 
than the ROPE pulse. 



optimal 

pseudopsectral 




FIG. 2: The efficiency of the transfer xi ~* Xi in system (|16p 
achieved by the pseudospectral method, as a function of the 
relaxation parameter ^ in the range [0, 1]. The theoretically 
calculated maximum efficiency given by (|17p is also shown. 



2. Spin Pair with Cross-Correlated Relaxation 

If DD-CSA cross-correlated relaxation cannot be ne- 
glected, the master equation as in (|15p is then modified 
to incorporate it as [10] 

p = - lJ[2Iizl2z,p] - kooWlzhz, Wlzl2z,p]] 

- khsAihz, [hz,p]] - kl;s^[l2z, [l2z,p]] 

- khn/csAi'^Iizhz, [hz,p]] 

DD/CSA [ 



[2/l./2.,[/2.,p]], 



where fc_D_D, ^^5^1, k'^sA ^'^^ auto- relaxation rates due to 
DD relaxation, CSA relaxation of spin Ii , CSA relaxation 
of spin I2 and k^j^^^g^, k'^j^^^g^ are cross-correlation 
rates of spins Ii and I2 due to interference effects between 
DD and CSA relaxation mechanisms. 
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(a) ROPE Control Amplitude (b) PS Control Amplitude 
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t t 



(c) ROPE Trajectories (d) PS Trajectories 

FIG. 3: Pseudospectral controls (panel b) and state trajec- 
tories (panel d) are compared to analytic ROPE controls 
(panel a) and trajectories (panel c) for ^ = 1. 
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FIG. 4: The efficiency of the transfer xi — > xe in system (|18p 
achieved by the pseudospectral method, as a function of the 
relaxation parameter in the range [0, 1], with = 0.75^0. 
The theoretically calculated maximum efficiency given by (|19p 
is also shown. 



Using this master equation, we can find the foUowing 
equations for the ensemble averages 
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where xi 



{2Ilj2z), 



(hy), Xi 



(18) 



^csa)/J^ 



(hz), X2 

X5 = {2Ii^l2z), X6 = i^hzhz), = {koD 

= ^DD/cSyi/"^ ^'^'^ ui{t),U2{t) are the available con- 
trols as before. Starting from x(0) — (1,0,0,0,0,0)^, 
we want to design ui(t) and U2(t) that maximize xe{T) 
subject to (fTS]) . 

This problem has also been solved analytically and the 
analytical pulse was denoted as CROP It is shown 
there that the maximum achievable value of a:6, i.e., the 
efficiency 772 of the transfer is given by the same formula 
as before 



but now 



(19) 



(20) 



Using the pseudospectral method introduced in this 
paper, we calculated numerically optimal controls 





0.5 1 1.5 2 2.5 3 



(a) CROP Control Amplitude (b) PS Control Amplitude 

FIG. 5: The CROP derived control pulse (panel a) is com- 
pared to the pseudospectral control pulse (panel b) for = 1 
and = 0.75. 



ui{t),U2{t) for various values of over [0,1] and with 
£_c = 0.75^a, to maximize the corresponding achievable 
values of X6{T). Using the same Matlab program and 
KNITRO solver, the optimal control problem was ap- 
proximated by 25 {N — 24) nodes with a free termi- 
nal time (maximum T = 5). A similar constant initial 
guess was used and the cost function tolerance was set to 
1 X 10^^. In Fig. [3]we plot the values of xq (T) achieved by 
the pseudospectral method and the maximum efficiency 
given in pop . Again, an excellent agreement is observed. 
The CROP and pseudospectral control pulse amplitudes 
are plotted in Fig. [51 
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B. Three Spin Chain 



The next open quantum system that we consider is a 
three spin chain (spins Ii, 12,13) with equal scalar cou- 
plings between nearest neighbors. In a suitably cho- 
sen (multiple) rotating frame, which rotates with each 
spin at its resonance (Larmor) frequency, the Hamilto- 
nian Hf that governs the free evolution of the system 
is Hf = V2J{Iizl2z + I2zl3z)- The common coupling 
constant is written in the form ^/2J for normalization 
reasons. As in the first example described in Section 
IIII A 11 we neglect cross-correlated relaxation and focus 
on slowly tumbling molecules in the spin diffusion limit. 
The corresponding master equation is 

p - iV2J[Iizl2z + hzhz,p] ~ kooWizhzj Wizhz^p]] 

- kooWizhz, ['^hzhz, p]] - kooi'^hzhz, ['^hzhzyp] 

- khsAlhz, [hz,p]] - klsA{l2z, [l2z,p]] 
-kl,sA[hz,[hz,p]]. (21) 

For this system we examine the polarization transfer from 
spin Ii to spin I3, Iiz I^z- This transfer is suitably 
done in three steps, hz ^ ^hzhz ~^ "^hzhz hz- 
The first and last steps are similar to the previously ex- 
amined spin transfer, thus we concentrate on the inter- 
mediate step 2/12/22 — 2/22/32, which corresponds to 
p(0) = 2/12/2Z and O — 2l2zl3z- Similarly, from the 
master equation (|2T|) . we can derive the associated dif- 
ferential equations which describe the time evolution of 
the expectation values of operators participating in this 
transfer. 



(22) 



where xi {^hzhz), ^2 = C^hzhx), = 

{V2{2hzl2yhz + l2y/2)), X4 = - {^hxhz) , X5 = 

{2l2zhz}, C = {2kDD + kcsA)/J UM- Transverse relax- 
ation rate is normalized with respect to the scalar cou- 
pling J and the normalized relaxation parameter is ^. 
The control function u{t) = LOy{t)/ J is the y-component 
of the applied magnetic field. The x-component creates 
an equivalent path from x\ to X5 and need not be con- 
sidered [13 ■ 

The corresponding optimal control problem is to find, 
starting from x{Q) = (1,0,0,0,0)-^, u{t) that maximizes 
X5{T) subject to (P^ . It has been shown in [l3| that a 
strict upper bound, 773, of the maximum achievable value 
of X5 is characterized by ^, 
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(23) 




This bound was used in ,17:] to evaluate the performance 
of a state-of-the-art gradient ascent algorithm for calcu- 
lating optimal u(t). 



FIG. 6: The efficiency for transfer xi X5 in system (|22p 
achieved by the pseudospectral method, as a function of the 
relaxation parameter ^ in the range [0, 1]. The correspond- 
ing numerical results achieved by a state-of-the-art gradient 
ascent algorithm as well as the theoretically calculated up- 
per bound for the maximum efficiency, given by (|23p . are also 
shown. 



In this paper, the optimal control u(t) was calculated 
by the pseudospectral method for various values of ^ 
in the range [0, 1] to maximize the corresponding val- 
ues X5{T), where the optimal control problem was ap- 
proximated by 25 {N — 24) nodes with a free termi- 
nal time (maximum T — 10). A similar constant initial 
guess was used and the cost function tolerance was set 
to 1 X 10~^. The values of Xr^{T) achieved by the pseu- 
dospectral method are displayed in Fig. [6l The corre- 
sponding numerical results by the state-of-the-art gradi- 
ent ascent algorithm used in 17] and the analytical effi- 
ciency upper bound given in (|23p are also shown in the 
same figure. Observe the excellent agreement between 
the efficiencies of the two numerical methods, lower of 
course than the analytical upper bound. In Fig. [7| we 
compare the gradient control pulse with the pseudospec- 
tral control pulse for the case ^ = 1. 



IV. CONCLUSION 

In this paper, we presented a method based on 
pseudospectral approximations to effectively discretize 
and solve optimal control problems associated with pulse 
sequence design for open quantum mechanical systems. 
Examples from NMR spectroscopy in liquid solutions 
evidenced the fiexibility and efficiency of the proposed 
methods. In these examples, pseudospectral methods 
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FIG. 7: The gradient derived control pulse (panel a) is com- 
pared to the pseudospectral control pulse (panel b) for ^ = 1. 

generated pulses that achieve performance similar to 
that of analytical methods, and it should be noted 
that these "approximated" optimal pulses found by 
pseudospectral methods are always smooth. A strength 
of pseudospectral methods is that they provide a robust 
technique which is easily extendible and implementable. 
In addition, it has been shown empirically that these 
methods have exponential convergence rates [33 |. 



while state-of-the-art algorithms like gradient methods 
typically evidence linear convergence. Pseudospectral 
methods provide a universal tool for solving pulse 
design problems for dissipative quantum systems coming 
from different physical contexts (NMR, MRI, Quantum 
Optics etc). Some immediate extensions of the methods 
presented here include considering the optimal pulse 
design problem of steering a family of open or closed 
quantum systems with different dynamics. A concrete 
example is to consider a family of coupled spin systems 
where each one is characterized by a different relaxation 
rate, e.g., ^ can take values from a positive interval 
[^1,^2]- In this case, the optimal control problem would 
seek to accomplish the maximum transfer, 77, over all 
systems, namely, to maximize J^^ viO^^ subject to the 
system evolution as in (|16p . (|18p . or ((22l) . Experimental 
verification of the performance of the pseudospectral 
pulses are also of keen interest and are currently being 
pursued. 
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